function [fm_tab] = firstmotion(traces,... traces:irisFetch trace struct
    T,... %periods
    starttime,...
    endtime,...
    minTseprat,...
    minSTALTA,...
    minampstd,...
    fm_LTA_window_cycles,... %number cycles for STA/LTA first motion pick
    fm_minampfrac,... %min amp for first motion pick relative to max (prevents picking artifacts)
    fm_tab) %optional
% calculate first motion of specified periods at each channel
% Josh Crozier 2019
% !!! currently only works when all traces have same sample freq !!!

fmplot = true;
VoicesPerOctave = 4;
WaveletParameters = [3,6];
%recommend 4-9 in order to resolve first peak, higher number increases
%larger artificial precursor oscilations
ExtendSignal = false;

%% main
if ~exist('fm_tab','var')
    fm_tab = table();
end

for Tind = 1:length(T)
    %window of interest
    local_starttime = max([starttime(Tind),traces(1).startTime]);
    startind = floor(seconds(local_starttime - traces(1).startTime)...
        *traces(1).sampleRate) + 1;
    local_endtime = min([endtime(Tind),traces(1).startTime...
        + seconds(traces(1).sampleCount/traces(1).sampleRate)]);
    endind = ceil(seconds(local_endtime - traces(1).startTime)...
        *traces(1).sampleRate) - 1;
    sampleCount = endind - startind + 1;
    if fmplot
        data_mat = zeros(length(traces),sampleCount);
    end
%     wt_mat = zeros(length(traces),sampleCount);
    iwt_mat = zeros(length(traces),sampleCount);
%     xcorr_mat = zeros(length(traces),sampleCount);
%     wvdcell = cell(length(traces));
    
% %     % bandpass filter
%     Fstophigh = 1/(minTseprat*seconds(T(Tind)))/2;     %Stopband Frequency
%     Fpasshigh = 1/(minTseprat*seconds(T(Tind)));     %Passband Frequency
%     Dpass  = 6;  % Passband Ripple
%     Dstop = 20;     %Stopband Attenuation
%     Fpasslow = minTseprat/seconds(T(Tind));     %Passband Frequency
%     Fstoplow = 2*minTseprat/seconds(T(Tind));     %Stopband Frequency
%     filt = designfilt('bandpassiir','Samplerate',traces(1).sampleRate,...
%         'StopbandAttenuation1',Dstop,...        
%         'StopbandFrequency1',Fstophigh,...
%         'PassbandFrequency1',Fpasshigh,...
%         'PassbandFrequency2',Fpasslow,...
%         'StopbandFrequency2',Fstoplow,...
%         'StopbandAttenuation2',Dstop,...
%         'PassbandRipple',Dpass,...
%         'MatchExactly', 'passband','DesignMethod','cheby2');
%     %fvtool(filt)

%     Dpass  = 4;  % Passband Ripple
%     Dstop = 20;     %Stopband Attenuation
%     Fpasslow = minTseprat/2/seconds(T(Tind));     %Passband Frequency
%     Fstoplow = minTseprat/seconds(T(Tind));     %Stopband Frequency
%     filt = designfilt('lowpassfir','Samplerate',traces(1).sampleRate,...
%         'StopbandAttenuation',Dstop,...  
%         'PassbandFrequency',Fpasslow,...
%         'StopbandFrequency',Fstoplow,...
%         'PassbandRipple',Dpass,...
%         'DesignMethod','kaiserwin',...
%         'ScalePassband',false);
%     filtdelay = round((filtord(filt)/2));

    if fmplot
        % bandpass filter
        Fstophigh = 1/(minTseprat*seconds(T(Tind)))/4;     %Stopband Frequency
        Fpasshigh = 1/(minTseprat*seconds(T(Tind)));     %Passband Frequency
        Dpass  = 6;  % Passband Ripple
        Dstop = 20;     %Stopband Attenuation
        Fpasslow = minTseprat/seconds(T(Tind));     %Passband Frequency
        Fstoplow = 2*minTseprat/seconds(T(Tind));     %Stopband Frequency
        filt = designfilt('bandpassfir','Samplerate',traces(1).sampleRate,...
            'StopbandAttenuation1',Dstop,...        
            'StopbandFrequency1',Fstophigh,...
            'PassbandFrequency1',Fpasshigh,...
            'PassbandFrequency2',Fpasslow,...
            'StopbandFrequency2',Fstoplow,...
            'StopbandAttenuation2',Dstop,...
            'PassbandRipple',Dpass,...
            'DesignMethod','kaiserwin',...
            'ScalePassband',false);
    %     filtdelay = round(mean(grpdelay(filt,linspace(Fpasshigh,Fpasslow,10),traces(1).sampleRate)));
        filtdelay = round((filtord(filt)/2));
    %     fvtool(filt)

    %     taperwin = tukeywin(sampleCount,0.1);
    end
    
    periodrange = [T(Tind),T(Tind)*2];
        
    for tracein = 1:length(traces)
%     try
%         column = strcat(traces(tracein).network,'_',...
%             traces(tracein).station,'_',...
%             traces(tracein).channel);

    if fmplot
            data_mat(tracein,:) = filter(filt,detrend(traces(tracein).data(startind:endind)));%.*taperwin);
            data_mat(tracein,:) = [data_mat(tracein,filtdelay+1:end),zeros(1,filtdelay)];
    %         data_mat(tracein,:) = smooth(traces(tracein).data(startind:endind),...
    %             round(seconds(T(Tind)*traces(tracein).sampleRate/2)), 'lowess');
    end
        
        [wt,Twt] = cwt(traces(tracein).data(startind:endind),'morse',seconds(1/traces(1).sampleRate),...
            'ExtendSignal',ExtendSignal,...
            'periodrange',periodrange,...
            'VoicesPerOctave',VoicesPerOctave,...
            'WaveletParameters',WaveletParameters);
%         [~,bestin] = min(abs(T(Tind)-Twt));
%         wt_mat(tracein,:) = wt(bestin,:);
        
        wt(2:end,:) = 0;
        iwt = icwt(wt,'morse',Twt,...
            [min(Twt),max(Twt)],...
            'WaveletParameters',WaveletParameters);
        iwt_mat(tracein,:) = iwt;

        %cross correlation wavelet
%         Ts = seconds(1/traces(tracein).sampleRate);
%         twave = Ts*(0:endind-startind);   
%         wave = sin(2*pi/seconds(T(Tind))*seconds(twave))...
%             .*exp(-twave*2*pi/(Q(Tind)*T(Tind)));
%         
%         %cross correlate
%         corrdata = traces(tracein).data(startind:endind);
%         corrstart = 0;%floor(seconds(T(Tind))/seconds(Ts));
%         xcorrvals = xcorr(corrdata,wave);
%         corrdata = xcorrvals(sampleCount-corrstart:end-corrstart);
%         xcorr_mat(tracein,:) = corrdata;
        
        %wigner ville
        %{
        x = traces(tracein).data(startind:endind);
        if mod(length(x),2) == 1
            x = x(1:end-1);
        end
        N = length(x);
        %t = 0:seconds(Ts):seconds(Ts)*N;
        xmirror = [fliplr(x'),x',fliplr(x')];
        wvd = zeros(size(x));
        mlist = -N:2:N;
        k = (seconds(Ts)*N)/seconds(T(Tind));
        for nind = 1:length(x)
            n = N + nind;
            temp = zeros(size(mlist));
            for mind = 1:length(mlist)
                m = mlist(mind);
                temp(mind) = xmirror(n + m/2)*xmirror(n - m/2)...
                    *exp(-1i*2*pi*k*m/N);
            end
            wvd(nind) = sum(temp);
        end
        %}
%         
%         [dwv,fwv,twv] = wvd(traces(tracein).data(startind:endind), 1/seconds(Ts));
%         dwv = abs(dwv);
%         [~,in] = min(abs(1./fwv-seconds(T(Tind))));
%         dwv = dwv(in,:);
%         wvdcell{tracein} = dwv;
        
    end %end loop over traces
    
    if fmplot
        %stack data
        data_stack = sum(abs(data_mat),1);
        total_std = std(data_mat(:));
        data_scaled_std = std((abs(data_mat)),1)./mean(abs(data_mat),1);
        data_atan_std = std(atan(abs(data_mat)*8*pi/total_std),1);
    end
    
    %stack wavelet
%     wt_stack = sum(abs(wt_mat),1);
% %     wt_angle_stack = sum(wrapToPi(angle(wt_mat)),1);
%     totalwt_std = std(abs(wt_mat(:)));
%     wt_scaled_std = std((abs(wt_mat)),1)./mean(abs(wt_mat),1);
%     wt_atan_std = std(atan(abs(wt_mat)*8*pi/totalwt_std),1);
    
    %stack inverse wavelet
    iwt_stack = sum(abs(iwt_mat),1);
%     totaliwt_std = std(abs(iwt_mat(:)));
%     iwt_scaled_std = std((abs(iwt_mat)),1)./mean(abs(iwt_mat),1);
%     iwt_atan_std = std(atan(abs(iwt_mat)*8*pi/totaliwt_std),1);
    
    %stack xcorr
%     xcorr_stack = sum(abs(xcorr_mat),1);
%     totalxcorr_std = std(xcorr_mat(:));
%     xcorr_scaled_std = std((abs(xcorr_mat)),1)./mean(abs(xcorr_mat),1);
%     xcorr_atan_std = std(atan(abs(xcorr_mat)*8*pi/totalxcorr_std),1);
    
%     %stack wvd
%     wvd_stack = wvdcell{1};
%     for tracein = 2:length(traces)
%         wvd_stack = wvd_stack + wvdcell{tracein};
%     end

    %% find onset
    %long term average and std
    LTA_window = ceil(fm_LTA_window_cycles*seconds(T(Tind))*traces(1).sampleRate);
%     LTA_mat = movmean(abs(data_mat),[LTA_window,0],2);
%     LT_std_mat = movstd(abs(data_mat),[LTA_window,0],0,2);
%     LTA_stack = sum(LTA_mat,1);
%     LT_std_stack = sum(LT_std_mat,1);
    LTA_iwt_mat = movmean(abs(iwt_mat),[LTA_window,0],2);
    LT_iwt_std_mat = movstd(abs(iwt_mat),[LTA_window,0],0,2);
    LTA_iwt_stack = sum(LTA_iwt_mat,1);
    LT_iwt_std_stack = sum(LT_iwt_std_mat,1);
    
    %criteria for possible onsets
    A_above_LTA = iwt_stack - LTA_iwt_stack;
    STALTA = iwt_stack./LTA_iwt_stack;
    iwt_range = range(iwt_stack(LTA_window+1:end));
    iwt_min = min(iwt_stack(LTA_window+1:end));
    maxin = islocalmax(A_above_LTA);
    maxin(maxin > LTA_window+1) = false; %trim start
    possiblein = maxin & ...
        A_above_LTA > minampstd*LT_iwt_std_stack & ...
        STALTA > minSTALTA & ...
        iwt_stack > iwt_range*fm_minampfrac + iwt_min...
            + min(iwt_stack(LTA_window+1:end));
    possiblein = find(possiblein);
    %pick first onset meeting all criteria
    if numel(possiblein) > 0
        onsetin = possiblein(1);
        %start time
        fm_tab{Tind,'fm_start'} = local_starttime + ...
            seconds((onsetin-1)/traces(1).sampleRate);
        fm_tab{Tind, 'fm_A_above_LTA'} = A_above_LTA(onsetin);
        fm_tab{Tind, 'fm_STALTA'} = STALTA(onsetin);
        fm_tab{Tind, 'fm_Afrac'} = iwt_stack(onsetin)/iwt_range;
        %find direction of each
        for tracein = 1:length(traces)
            column = strcat(traces(tracein).network,'_',...
                traces(tracein).station,'_',...
                traces(tracein).channel);
            fm_tab{Tind, strcat('fm_',column)} = iwt_mat(tracein,onsetin);
        end
    else
        onsetin = 1;
    end
    
    %% debug plot
    if fmplot
        figure(16),clf
        plottime = linspace(local_starttime,local_endtime,sampleCount);

        subplot(2,1,1)
%         plot(plottime,abs(data_mat(1,:))), hold on
%         for tracein = 2:length(traces)
%             plot(plottime,abs(data_mat(tracein,:)))
%         end
        plot(plottime,data_stack,'k'), hold on
        %title(strcat('T=',num2str(seconds(T(Tind)))));
        for ii = 1:length(possiblein)
            xline(plottime(possiblein(ii)),'m');
        end
        xline(plottime(onsetin),'c');
        title(strcat('Bandpass:',num2str(1/Fpasslow,'%.1f'),...
            '-',num2str(1/Fpasshigh,'%.1f'),' s'))
        ylabel('Stacked Amplitude')

    %     subplot(6,1,2)
    % %     plot(plottime,abs(LTA_mat(1,:))), hold on
    % %     for tracein = 2:length(traces)
    % %         plot(plottime,abs(LTA_mat(tracein,:)))
    % %     end
    % %     plot(plottime,LTA_stack,'k')
    %     plot(plottime, data_stack-LTA_stack,'k')
    %     title('LTA');
    %     for ii = 1:length(possiblein)
    %         xline(plottime(possiblein(ii)),'m');
    %     end
    %     xline(plottime(onsetin),'c--');

    %     subplot(6,1,3)
    %     plot(plottime,data_atan_std,'k')
    %     hold on
    %     plot(plottime,data_scaled_std,'b')
    %     title('std and std of arctan')
    %     legend('atan','std')
    %     for ii = 1:length(possiblein)
    %         xline(plottime(possiblein(ii)),'m');
    %     end
    %     xline(plottime(onsetin),'c--');

    %     subplot(6,1,3)
    %     plot(plottime,abs(wt_mat(1,:))), hold on
    %     for tracein = 2:length(traces)
    %         plot(plottime,abs(wt_mat(tracein,:)))
    %     end
    %     plot(plottime,wt_stack,'k')
    %     title('wt');  
    %     
    %     subplot(6,1,4)
    %     plot(plottime,wt_atan_std,'k')
    %     hold on
    %     plot(plottime,wt_scaled_std,'b')
    %     title('std and std of arctan')
    %     legend('atan','std')

        subplot(2,1,2)
%         plot(plottime,abs(iwt_mat(1,:))), hold on
%         for tracein = 2:length(traces)
%             plot(plottime,abs(iwt_mat(tracein,:)))
%         end
        plot(plottime,iwt_stack,'k'), hold on
        title(strcat('Wavelet Filter',{' '},'T=',num2str(seconds(T(Tind))),' s'));  
        for ii = 1:length(possiblein)
            xline(plottime(possiblein(ii)),'m');
        end
        xline(plottime(onsetin),'c');
        ylabel('Stacked Amplitude')

    %     subplot(2,1,2)
    % %     plot(plottime,abs(LTA_iwt_mat(1,:))), hold on
    % %     for tracein = 2:length(traces)
    % %         plot(plottime,abs(LTA_iwt_mat(tracein,:)))
    % %     end
    % %     plot(plottime,iwt_stack,'k')
    %     plot(plottime, iwt_stack-LTA_iwt_stack,'k')
    %     title('iwt LTA');  
    %     for ii = 1:length(possiblein)
    %         xline(plottime(possiblein(ii)),'m');
    %     end
    %     xline(plottime(onsetin),'c--');

    %     subplot(6,1,6)
    %     plot(plottime,iwt_atan_std,'k')
    %     hold on
    %     plot(plottime,iwt_scaled_std,'b')
    %     title('std and std of arctan')
    %     legend('atan','std')
    %     for ii = 1:length(possiblein)
    %         xline(plottime(possiblein(ii)),'m');
    %     end
    %     xline(plottime(onsetin),'c--');

    %     subplot(6,1,5)
    %     plot(plottime,abs(xcorr_mat(1,:))), hold on
    %     for tracein = 2:length(traces)
    %         plot(plottime,abs(xcorr_mat(tracein,:)))
    %     end
    %     plot(plottime,xcorr_stack,'k')
    %     title('xcorr');  
    %     
    %     subplot(6,1,6)
    %     plot(plottime,xcorr_atan_std,'k')
    %     hold on
    %     plot(plottime,xcorr_scaled_std,'b')
    %     title('std and std of arctan')
    %     legend('atan','std')

    %     subplot(5,1,5)
    %     plot(local_starttime + seconds(twv),wvdcell{1}), hold on
    %     for tracein = 2:length(traces)
    %         plot(local_starttime + seconds(twv),wvdcell{tracein}), hold on
    %     end
    %     plot(local_starttime + seconds(twv),wvd_stack,'k')
    %     title('wvd');
        %}
    end
               
%     catch ME
%         disp(ME.message)
%         traces(tracein)
%         tracein
%         T
%         any(~isfinite(traces(tracein).data))
%     end
end %end loop over traces

end

